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MULTI-COMPONENT SEISMIC DATA PROCESSING 

The present invention generally relates to apparatus and 
methods, for processing multi- component (3C/4C) seismic data 
5 based on the use of receiver functions, more specifically 

« 

two-dimensional receiver functions-. It also pertains to the 
use of receiver functions to process and interpret seismic 
signals to derive an earth image away from the near- surface 
structure at the receiver location. 

10 

BACKGROUND OF THE INVENTION 

t 

* \ 

The processing and interpretation of multicomponent -(3C/4C) 

w 

seismic data, acquired directly at the seafloor or during 3C 
land seismic surveys is compromised by the effects that the 

15 shallow subsurface has on the deeper reflected wavefield. 
The near surface is generally associated with low, laterally 
varying shear wave velocities and on land, the P-wave 
velocity can also be low. These properties often lead to 
• large P- and S-wave traveltime perturbations in the deeper 

20 reflected wavefield which vary from receiver to receiver. In 

■ 

addition, there are suggestions that ' scattering and 
anelastic attenuation (especially shear) impacts amplitudes 
and recorded wavef orms as well . 

25 The receiver function methodology has it roots in earthquake 
seismology, where it was developed to investigate the 
structure of the crust and upper mantle using multicomponent 
teleseismic body wave data. Certain aspects of receiver 
functions are described for example by T. Ryberg and M. 

30 Weber in Geophys. J. Int. 141 (2000), 1-11. Reference to 
aspects of receiver function can also be found in the 
published patent applications GB 2384557 and WO 02/059647. 
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« 

While trying to image discontinuities in upper mantle using 
receiver functions, it was realised that the change in the 
travel time difference between the P-wave and the PS- 
converted wave, as a function of the angle of incidence, 
5 could no longer be neglected- T. Ryberg and M. Weber 
proposed to calculate a velocity spectrum stack (VSS) based 
on the change in the traveltime difference (hereafter called 
moveout) in a single layer to find optimal velocities for 

stacking of receiver functions. Although this approach gives 

♦ 

10 good results when the structure is relatively homogeneous 
above the converting interface, it lacks a theoretical basis 
when applied to layered media with significant other 
discontinuities above the converting interface, i.e. when 
ray bending is likely to produce significant deviations from 

15 the traveltime difference equation for a single layer. This 
was partly recognised and, subsequently, it was suggested to 
substitute the average vertical P- and S-wave slowness for a 
stack of layers, into the equation for a single layer. It 
was also proposed to circumvent this problem by assuming a 

20 reference model (providing an initial vp(z) and vs(z)) and a 

• ■ 

corresponding one-parameter family of related models by 
multiplying the reference model by a fraction close to one, 
and obtain the difference in traveltime as a function of 
slowness by ray-tracing through the perturbed reference 
25 model. This approach severely restricts the number of free 
parameters . 

On the other hand, normal moveout correction (NMO) , velocity 
analysis and stacking of multichannel data has a long 
30 history in exploration and production seismic and the theory 
is well developed for horizontally layered media. It was 
demonstrated how to derive a power series approximation of 
the squared reflection traveltime as a function of even 
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powers of offset (i.e. of the form: t 2 «Ci+c 2 -x 2 +c 3 -x 4 +..., where 
t denotes two-way traveltime, x denotes offset and the 
coefficients depend on layer thicknesses and velocities of 
the medium) . Hyperbolic approximations were used, truncating 
the infinite series after the second term and defined the 
rms-velocity as the square-root of the inverse of 
coefficient c 2 . It is also known that the slope of the x 2 -t 2 
curve at x=0 yields the inverse of the squared rms-velocity 
and how from the rms-velocity at two consecutive depth 
levels the interval velocity between them can be calculated. 
Taner and Koehler in: Geophysics 34 (1969), 859-881 also 
introduced the velocity spectrum stacking technique using a 
multichannel coherence measure called semblance. The work by 
Taner and Koehler was generalised for PS-converted waves by 
Tessmer and Behle in Geophys . Prosp. 36 (1988),. 671-688, who 
also derive a Dix-Krey type formula, relating rms-velocities 
for PS-converted waves to products of P- and S-wave interval 
velocities in each layer. 

In view of the above state of art, it is an object of the 
present invention to extend and improve the use of receiver 
functions to process and interpret seismic data to derive 
images of an earth. It is further object of the invention to 
^ determine velocities or velocity models from such receiver 
functions . 

SUMMARY OF THE INVENTION 

This invention describes systematic approximations of the 
difference in traveltime between a seismic signal in an 
unchanged mode and a corresponding mode-converted signal, 
such as P-waves (or S-waves) and PS-waves (or SP- waves) , 
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converted upon reflection or transmission at a particular 
interface in a horizontally layered medium or a multi- 
layered earth model. As such, the relations of the present 
invention can be used to approximate the moveout of events 
in receiver functions, which correspond to traveltime 
differences between non-converted and converted waves. The 
receiver function can thus be corrected for moveout in the 
slowness domain or in the time- space domain. The moveout 
correction overcomes limitations of known attempts to use 
receiver function to interpret an earth, which attempts were 
limited to shallow subsurface layers or cases of low- 
velocity layers (or vertical or near-vertical raypaths) . 

The receiver functions are well known seismology but rarely 
used in seismic operations related to the exploration and 
monitoring of hydrocarbon reservoirs. The receiver function 
are broadly defined as the result of a deconvolution or 
crosscorrelation between mode- converted events at a receiver 
location, thereby representing an ideally source- independent 
response of the earth layers beneath the receiver or more 
generally beneath the location where the deconvolution or 
crosscorrelation is performed. 

Receiver function can be derived in a variety of ways mainly 
distinguished by the use of different components of the 
wavefield or a different set of coordinates. They include: 

- the deconvolution or crosscorrelation of vertical and 
horizontal components (both in-line and cross-line) of 
particle displacement / velocity / acceleration; 

- the deconvolution or crosscorrelation of vertical and 
radial /transverse components (or any component resulting 
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from rotation of linearly independent horizontal 
measurements into a particular coordinate system) ; 

- the deconvolution or crosscorrelation in the LQT system 
often applied in seismology with the L component (pointing 

5 in the direction of propagation of an incoming P-wave, below 
the free-surface) and the Q- and T-components (mutually 
orthogonal and orthogonal to the L component, hence pointing 
in the direction of particle motion of the SV and SH-waves 
respectively) ; 

10 - the deconvolution or crosscorrelation of any components 
resulting from rotation of 3 mutually orthogonal 
measurements of particle motion in a particular coordinate 
system, with the aim of better separating P- and S-waves 

■ 

before calculating the receiver function; 

15 - the deconvolution or crosscorrelation of any P- and S-wave 
quantities, either measured directly (e.g. using special 
receiver configurations measuring divergence or curl, or 
special sensors (e.g. a pressure sensor etc.) or resulting 
from wavefield decomposition of mutually orthogonal 

2 0 measurements ; 

- the deconvolution or crosscorrelation of complete 
recordings as well as selected parts of recordings; 

- the equivalents of deconvolution and crosscorrelation in 
the frequency (or frequency-horizontal wavenumber or 

25 frequency- slowness) domain; 

- 2D- or 3D- deconvolution / crosscorrelation involving 1- 
or 2- spatial axes in the deconvolution / correlation 
process; or 

- deconvolution or crosscorrelation of P- and PS-converted 
30. waves as well as S- and SP-converted waves (i.e. all 

combinations of converted and unconverted waves) 
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The receiver function deconvolution/crosscorrelation is not 
limited to a physical surface: it may be done everywhere in 

» 

the medium (e.g. as part of a wavefield extrapolation 
prpcedure) . 

5 

The receiver function deconvolution/crosscorrelation 
procedure of the present invention is not necessarily 
limited. to taking pairs of traces from the same receiver, it 
can be the result of a combination of particular traces 
10 containing P- and S-wave events from different receiver 
locations . 

In a, first aspect the traveltime difference is expanded into 
increasing even powers of slowness, leading to a generalised 
15 small-slowness approximation of receiver function moveout in 
the time-slowness domain. This result is more genera.1 than 
previous results by Ryberg and Weber (2000) since it applies 
to a stack of layers and contains their results as the 
special case of a single layer over a half space. A pseudo 

2 0 rms-velocity is defined and this leads naturally to a Dix- 

Krey type relation for receiver functions, relating the 
pseudo rms-velocities for two consecutive interfaces to the 

a 

■ 

product of the P- and S-wave interval velocity between them. 

25 In a second aspect, a related expansion of the squared, 
traveltime difference is made into increasing even powers of 
tlie difference in horizontal travel distance between the P- 
and PS-converted wave (e.g., dt 2 =ci+ c 2 • dx 2 +c 3 • dx 4 +..., where dt 
denotes the difference in traveltime between the PS- 

3 0 converted wave and the P-wave from a particular interface, 

whereas dx denotes the difference in offset of emergence 
between the P- and PS-wave) . In connection with this 
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expansion, a new type of receiver function is also 
introduced, formed by two-dimensional deconvolution of (x,t) 
domain vertical component data (i.e. P-waves) out of the 
corresponding (x,t) radial component data ( PS - waves ) . This 

■ 

5 2D-receiver function, because of the added deconyolution in 
the spatial direction, gives precisely the time-lag between 
the PS-wave and the P-wave as a function of the distance- lag 
between them. The second (truncated) expansion hence can be 

■ 

used to NMO correct the new type of receiver function. 
10 Synthetic data generated with a reflectivity code is used to 
illustrate the moveout approximations . and the new type of 
receiver function . 

The derived series approximations (and corresponding two- 

15 term truncations) of travel time differences between PS- 
converted waves and P-waves from interfaces in horizontally 
layered media, being either a function of slowness or a 
function of the difference in horizontal distance traveled 
by the two phases, enable velocity analysis and subsequent 

20 normal moveout correction of events in receiver functions. 
The derived Dix-Krey type formulas for receiver functions 
relate the (pseudo) rm's -velocities for two consecutive 
interfaces in a layered earth to the product of P- and S- 
wave interval velocities between them. A new kind of 

25 receiver function is proposed, calculated by two-dimensional 
deconvolution of (x,t) domain vertical component data out of 
(x,t) radial component data, which gives information about 
the travel time difference between P- and PS-converted phases 
as a function of the difference in horizontal distance 

30 traveled by those phases. The second moveout approximation 
applies to events in this new type of receiver function. 
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Aspects of the invention further include a method of 
calculating the moveout (or change as a function of 
slowness) of the traveltime difference between a PS- 
converted wave and a P-wave, reflected/ transmitted from the 
same interface in a layered earth, by expanding the 
traveltime difference into an infinite series of increasing, 
even, powers of slowness. This approximation can be made 
increasingly accurate by including more terms in the 
expansion. Coefficients of higher-order terms can' be 
obtained by regression analysis of picked traveltime 
differences in the (T,p) domain receiver functions. 

Another aspect of the invention is a two-term truncation of 
the series expansion of the traveltime difference between a 
PS-converted wave and a P-wave reflected/transmitted from 
the same interface and the association of the coefficient 
multiplying the second term with a pseudo rms-velocity . This 
approximation expresses the moveout of the traveltime 
differences a in terms of a single profile or parameter of 
scalar quantities with depth rather than two profiles such 
as the P- and S-wave velocity profiles Vp(z) and Vs(z) . 

* 

Yet another aspect of the invention is a method of velocity 
analysis of receiver functions calculated by deconvolution 
or crosscorrelation of vertical component data (or any P- 
wave quantity) with radial component data (or any S-wave 
quantity) in the (x,p) domain, which uses the above two-term 
approximation and any multichannel coherency measure, 
similar to the velocity analysis method proposed by Taner 
and Koehler (1969) for P-wave data or other known methods. 

Yet another aspect of the invention is a method of 
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processing the (pseudo-) velocities produced by the method 
of receiver .function velocity analysis, as described above, 
to obtain the product , of p- and S-wave propagation 
velocities in each layer of the medium, using the Dix-Krey 
type relations developed for the (x,p) domain receiver 
functions . 

Yet another aspect of the invention is method of correcting 
the events in a (x,p) domain receiver function for the NMO 
of the travel time difference using the product of the P- and 
S-wave interval velocities obtained using the method 
described above and a priory knowledge of either the P- or 
S-wave velocity or velocity ratio in each layer and 
theoretical expressions or ray- tracing. 

Another aspect of the invention is a method of calculating 
two-dimensional (2D) receiver functions, by 2D (stabilised) 
deconvolution or crosscorrelation of (x,t) domain vertical 
component data (or any P-wave quantity) with the 
corresponding (x,t) domain radial component data (or any S- 
wave quantity) 

Yet another aspect of the invention is a method of 
calculating the moveout (or change as a function of the 
difference in horizontal • distance traveled) of the 
travel time difference between a PS-converted wave and a P- 
wave reflected/ transmitted from the same interface, by 
expanding the squared traveltime difference into an infinite 
series of increasing powers of the difference in horizontal 
distance traveled by both phases. Since the result is an 
asymptotic expansion (around dx=0) , including more terms in 
the expansion does not necessarily improve the accuracy. 
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The expansion is truncated one term before the smallest 
which guarantees that the error in the expansion is of the 
order of the first neglected term. 

Yet another aspect of the invention is a method of 

4 

calculating a similar asymptotic expansion as above around a 
different expansion point (e.g. dx=250m) while maintaining 
the possibility to interpret the first few coefficients as 
the vertical incidence traveltime and the rms-velocity 
respectively. 

Yet. another aspect of the invention is a two-term truncation 
of the series expansion of the squared traveltime difference 
as a function of increasing even powers of the difference in 
horizontal distance traveled by the two- phases and the 
association of the coefficient multiplying the second term 
with a rms-velocity. 

« 

Yet another aspect of the invention is a method of velocity 
analysis of the 2D-receiver functions formed by two- 
dimensional deconvolution using the two- term truncation of 
the traveltime difference as a function of the difference in 
horizontal distance traveled and a coherency measure, 
similar to the velocity analysis method proposed by Taner 
and Koehler (1969) for P-wave data. 

Yet another aspect of the invention is a method of 
processing the velocities produced by the method of receiver 
function velocity analysis to obtain the product of P- and 
S-wave propagation velocities in each layer of the medium, 
using the Dix-Krey type relations developed for the (x, t) 
domain two-dimensional receiver functions. 
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m 

Yet another aspect of the invention' is a method of 
correcting the events in a two-dimensional receiver function 
for the NMO of the travel time difference using the product 
5 of the P- and S-wave interval velocities obtained using the 
above-described method and a priori knowledge of either the 
P- or S-wave velocity or velocity ratio and theoretical 
expressions or ray- tracing. 

10 The methods of the invention, including the moveout 
corrected receiver functions or the velocities derived from 
the related Dix-Krey relations, can be applied to many known 
seismic data processing methods, including velocity 
analysis, moveout correction (NMO and/or DMO) , stacking and 

■ 

15 pre- or post-stack migration. 

These and other aspects of the invention will be apparent 
from the following detailed . description of non- limitative 
examples and drawings. 

20 

BRIEF DESCRIPTION OF THE DRAWINGS 

V 

FIG 1 is a schematic illustration of the phases and raypaths 
involved in the seismological receiver function 
setting. A plane P-wave is incident on a stack of n 

25 horizontal layers from below and mode-converts to 

Shear at discontinuities within, and at the base of 
the stack. P-waves are denoted by solid lines, S- 
waves by dashed lines. The data from each 
multi component station is processed separately and 

3 0 the extra time Th is treated implicitly in the 

traveltime difference calculation. 
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FIG 2 is a schematic illustration of the phases and 
raypaths involved in the reflection seismic receiver 
function setting. A P-wave is incident on a stack of 
n horizontal layers from above and mode-converts to 
5 Shear upon reflection at discontinuities within and 

at the base of the stack. P-waves are denoted by 
solid lines, S-waves by dashed lines. The data is 
recorded on an array of equidistant multicomponent 
receivers and the difference dx in horizontal 
10 distance graveled by the P- and PS -converted wave . 

from the same reflection/ conversion point can be 
treated explicitly, analogous to the difference in 
traveltime between the two phases . 

15 FIG. 3 shows pre-processed reflectivity data for the six- 

layer over a halfspace model (medium properties in 
Table 1 below) . The primary P- and PS-converted waves 
have been identified by ray-tracing through the model 
and are indicated in blue and red respectively. This 

20 is the data input for receiver function calculation 

in the slowness domain. 

FIG. 4A are receiver functions obtained by stabilised 
deconvolution of the upgoing P-waves out of the 

25 upgoing PS-converted waves shown in FIG. 3. The 

theoretical traveltime differences between the PS- 
converted waves and primary P-ref lections are 
calculated by ray- tracing and are shown in blue. The 
two- term approximation (eq. 5) using exact pseudo 

30 rms-velocities is shown in cyan. The black box 

denotes the area which is shown in detail in FIG. 5A. 
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FIG. 4B are NMO corrected receiver functions using the two- 
term approximation (eq. 5) and the linearly 
interpolated exact pseudo rms-velocity model . The 
predicted residual moveout is shown in blue. The 
black box denotes the area which is shown in detail 

* 

in FIG. 5B. 

■ 

FIG. 5A is a zoom- in on the area denoted by the black box in 
FIG 4A. The two main events correspond to the 
relative moveout between the P-reflection and its PS- 
conversion at the fifth and sixth interface. It is 
showing receiver functions obtained by stabilised* 
deconvolution of .the upgoing P-waves out of the 
upgoing PS-converted waves shown in FIG. 3. The 
theoretical traveltime differences between the PS- 
converted waves and primary P-ref lections are 
calculated by ray- tracing and are shown in blue. The 
two-term approximation (eq. 5) using exact pseudo 
rms-velocities is shown in cyan. 

* 

FIG'. 5B is a zoom- in on the area denoted by the black box in 
. FIG . 4B . The two main events correspond to the 
relative moveout between the P-ref lection and its 
PS-conversion at the fifth and sixth interface. It 
is showing NMO corrected receiver functions/ using 
the two- term approximation and the linearly 
interpolated .exact pseudo rms-velocity model. The 
predicted residual moveout is shown in blue. 
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FIG. 6A is a receiver functions calculated by two- 
dimensional stabilised deconvolution of (x,t) domain 
vertical component data out of the corresponding 
radial component data. This 2D-Receiver Function 
5 gives the traveltime difference between a P-wave 

reflection and its PS-conversion as a function of the 
difference in the horizontal travel distance'. The 
theoretical traveltime differences, calculated by 
raytracing, are indicated in blue. The two-term 
10 approximation, calculated from eq. 20 and the 

coefficients given in 25a, b using exact medium 
properties, is indicated in cyan. 



FIG. 6B shows again theoretical traveltime differences and 
15 the two- term approximation are indicated in blue and 

cyan respectively. Successive higher order 
approximations are shown in green. 



FIG. 7A is a zoom- in on the 2D-Receiver Functions and 
20 moveout approximations shown in FIG. 6A. The two main 

events give information about the difference in 
traveltime between the P-reflection and its PS- 
conversion at the fifth and' sixth interface as a 
function of the difference in horizontal distance 
25 traveled. Theoretical traveltime differences and the 

two- term approximation (eq.20) are indicated in blue 
and cyan respectively. 



FIG. 7B is a zoom- in on the 2D-Receiver Functions and 
30 moveout approximations shown in FIG. 6B. Successive 

higher order approximations are also shown,- in green 
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DETAILED DESCRIPTION AND EXAMPLES 

5 Series expansion of the traveltime difference as a function 
of slowness 

It is well known from global seismology that the difference 
in traveltime between a P-wave incident on a stack of n 
10 horizontal layers from below and the corresponding PS-wave, 
converted upon transmission into the stack, can be written: 

* = E K {jvs; 2 -p 2 - Jv P ? - P 2 ), (i) 

i 

15 where h k ,Vs k and Vp k denote the thickness, shear and 
compressional wave velocity of layer k respectively, and p 
the slowness of the incident plane-wave. Note that this 
formula implicitly takes into account the extra time (Th) it 
takes the incident plane-wave to reach the horizontally 

20 offset PS-conversion at the base of the stack. Hence this 
formula can be used to compare the arrival time of plane P- 
and PS-converted waves with the same slowness, arriving at 
the same receiver, as is normally done in seismology (see 
FIG 1) . Equation (1) contains terms of the form: (l-p 2 -^)* 6 , 

25 where v can denote either the P- or S-wave velocity. Such 
terms can be expanded into a Taylor series as follows: 

a-pV)* = i>,(/>vy = i-i(pV)-i(A 2 ) 2 -£(pV> 3 o 

* 

30 where the coefficients qj are given by: 

15 
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(3) 



Substituting equation (2) into equation (1), the traveltime 
5 difference can be written: 



00 n 



dt --SffiZw^ -wv ■ (4) 

7=0 



10 Equation (4) ( is an expansion of the traveltime diffe^nce 
between plane P- and PS- converted waves through an n- layered 
medium in increasing even powers of slowness. Note that the 
accuracy of the expansion is directly related to the 
accuracy of the substituted Taylor series approximation 

i5 (equation 2) , which means that the product p 2 v* should be 
small (at least smaller than 1, which means that the waves 
should be • nowhere evanescent) throughout the stack. 
Therefore the product of the highest P-velocity and the 
slowness determine the accuracy of equation (4) . If it is 

20 assumed that the slowness is small throughout the stack, we 
can truncate the infinite series given by eq. (4) after the 
second term and neglect terms of fourth order in slowness- 
and higher. This gives the small- slowness, or short-spread 
approximation : 

25 

« • 
dt = XhtW? -Vp^ + ^KiVpt -Vs k )p 2 . (5) 

*=1 *=1 



The first term in equation (5) is simply the difference in 
30 traveltime between a vertically incident P- and S-wave (p=0) 
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through the stack of layers. The second term, multiplying 
p , is not simply interpretable since it contains products 
of layer thicknesses and differences in velocities and has 
units [m /s] . However, maintaining the analogy with NMO 
corrections developed in exploration and production seismic, 
we call this the pseudo rms-velocity. Thus the traveltime 
difference at vertical incidence dt 0 and the pseudo rms- 
velocity v^ respectively can be defined as: 

r 

rf'o = E w; 1 - vp; 1 ), (6a) 

n 

Vnns = l2X <VPk "VS.). (6b) 

Note that for a single layer, the short-spread approximation 
given in equation (5) reduces to the equation given by 
Ryberg and Weber (2000). Hence, their result is a special 
case of the small-slowness relation proposed here. Note also 
that although equation (1) originates from earthquake 

seismology and is usually only applied in a* transmission 

« 

setting, it is also valid in a reflection setting, where 
plane P- and PS-waves aire considered, reflected and mode- 
converted at the same interface in a horizontally layered 
medium and recorded at the same receiver. 

A Dix-Krey relation for receiver function pseudo rms- 
velocities 

m 

The definition of a pseudo rms-velocity in equation (6b) 
allows the derivation of a Dix-Krey type formula that 
relates the pseudo rms-velocities to the product of P- and 
S-wave interval velocities, as will be shown below. A 
similar analysis as Tessmer and Behle (1988) is applied and 
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the thickness h k of layer k is expressed in terms of the 
one-way vertical travel times ' r p k and r s * and interval 
velocities in that layer: 



k=Wp&+ w*<>. (7) 

Also, the ratio of vertical one-way. travel times equals the 
inverse of the ratio of interval velocities: 



< v Pk 



(8) 



If equation (7) is subsituted . into the. definition of the 
pseudo rms-velocity (eq. 6b) and the identity in eq. (8) is 
used, one can derive the following expression for the pseudo 
rms-velocity at layers n and n-lx 



*L=VH?l-*D-Vs k Vp k , (9a) 



*=1 
n-l 



-VWk. * (9b) 

k=l ■ 



Subtracting equation (9b) from equation (9a) , one arrives 
at : 



Vs n Vp n =2 ^ — 



(10) 



Finally, using that the difference in the vertical one-way 
S- and P-wave traveltime through layer n is equal to the 
difference in vertical incidence traveltime differences 
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between layer n and layer n-1 



{rl-T p k ) = dt n Q -d%-\ (11) 



one arrives at a Dix-Krey type of relation that relates the 
pseudo rms-velocities to the product of P- and S-wave 
interval velocities : 



Vs n Vp n = 2 11113 ^ 



dtl - dt 0 



n-1 ■ 



(12) 



Equation (12) says that the product of P- and S-wave 
interval velocities in layer n is twice the ratio of the 
differences in pseudo rms-velocities and vertical incidence 
traveltime . differences . Equation (12) can be verified by 
directly substituting the definitions of the pseudo rms- 
velocities and vertical incidence traveltime differences 
(i.e. eq. 6a, b) . 



Expansion of the traveltime difference as a function of the 
difference in horizontal distance traveled 



In the previous sections, the horizontal distance between 
the point of P-wave transmission and PS-conversion at the 
base of the stack, but emerging at the same surface 
location, was taken into account implicitly in the 
derivation of the traveltime difference (equation 1) . The 
main (historical) reason for this is that it allows accurate 
comparison of P- and PS-wave traveltimes at a single 
multicomponent station (see FIG 1) . However, when the 
multi component data are recorded on a densely spaced array 
of receivers, such as in exploration and production seismic, 
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the P- and PS-wave energy originating from the same location 
on the converting interface are both recorded, although at 
horizontally offset locations in the receiver array- This 
allows us to consider the difference -in traveltime and in 
5 horizontal travel distance explicitly for an array of 
receivers {compare FIG.l and FIG. 2) 

Using elementary trigonometric relations and Snell's law, 
the explicit difference in traveltime between a plane P- and 
10 PS-wave of slowness p, transmitted and converted at the same 
location at the base of layer n, can be written (where h± 
,Vs± and Vp± denote the thickness, shear and compressional 
wave velocity of layer i respectively) : 



15 dt = Y t (Ts k -7}> k )='Z 



K 



K 



Vs k J\-p 2 Vsl Vp k j\-p 2 Vp 



2 
k 



(13) 



Similarly, the difference in horizontal distance traveled 
can be written: 



n 



n 



20 dx = Y J (Xp k -Xs k )=5r 



Jfc=l 



KVp* 



(14) 



25 



Now, following the analysis of Tessmer and Behle (1988) , 
both dt and dx are expanded into infinite series of even 
powers of p. Hereby the Taylor series expansion of the 
function (l-p 2 ^ 2 ) ~ H is used, where v can denote both the P- 
and S -wave velocity. 



(1-pV)-* = t^CpV)'- 1 =i + I(A 2 ) + ^|(pV) 2 + ^|( P V) 3 + ..., (15) 
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where the coefficients Qj are given by 



_ lx3x--x(27-3) 

q x =i, Qj-- — : — — — . (16) 

7 2x4x-x(27-2) 



Note that this expansion is the inverse of the expansion 
used in the previous sections (see equation 2) . Thus 
substituting equation (15) into equations (13) and (14) for 
10 both terms containing P- and S-velocities to obtain infinite 
series for dt and dx results in: 



n 



* - E<Z;2X<^" 3 -Vp 2 k J - 3 )(p 2J - 2 ) , (17a) 



n 



** = 2>,Z W*r x -Vp^XP 27 " 1 ). (17b) 

j=\ k=l 



To simplify the appearance of equations (17a, b) and 
subsequent derivations, the' following coefficients are 
defined: 



n 



a m =HK(V4 m ~ 3 -Vpl m - 3 ), (18a) 

20 b m =-q m a m+l , (18b) 

(18c) 



m "am m+l ' 
^ = 1m a m ■ 



Using these, equations (17a, b) become 



2*-2 

*=1 



OB 



(19a) 



(19b) 
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Equations '19a, b have the same form as Taner and Koehler 
(1969) derived for the traveltime and offset of a P-wave 
reflection in a horizontally layered medium and as Tessmer 
5 and Behle (1988) found for the traveltime and offset of PS- 
converted waves. This suggests that it is possible to apply 
the same methodology as they have done to find an expansion 
of the square of the traveltime difference into increasing 
even powers of the difference in horizontal traveldistance . 
10 Such an expansion would be useful if a new kind of receiver 
function could be obtained, where PS- converted waves are 
shifted both in time and space relative to the P-wave, 
measuring the traveltime difference as a function of the 

i 

difference in horizontal distance traveled by both phases. 

15 But this new receiver function is exactly the result of 
applying a two-dimensional (2D) deconvolution of the (x,t) 
domain vertical component data from the - (x,t) domain 
horizontal component data. The 2D or 3D deconvolution using 
one or two spatial axes yields directly the distance lags dx 

20 in the same way as the deconvolution over time yields the 
time-lags dt. Thus, the objective is to express the 
traveltime difference dt as: 



The explicit form of equation (20) can be found by squaring 
the power series for dt (equation 19a) , successive higher 
even powers of dx (equation 19b) and ordering terms of equal 
30 powers in slowness p. Starting by squaring dx: 



dt 2 =c x +c 2 dx 2 +c 3 dx 4 +<? 4 -dx 6 + 



(20) 
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dx 1 = 



pZhp 2 *- 2 = P % +b 2P 2 +b 3P * +b iP 6 +-) 2 = 



(21) 



p 2 (^ 2 +(bA +b 2 b l ) P 2 +(b l b 3 +b 2 +bj> l ) P * +(*,& 4 +b 2 b 3 +b 3 b 2 +b A b0p 6 +~)=P 2y ZB klP 2h - 2 , 



*=i 



where , 



5 B kl = + fcA-i + • • • + h-A + Kbx • 



(22) 



10 



The higher, even powers of this series can be found by 
recursively applying equation (22) to find the coefficients. 
This leads to the following result: 



(23) 



*=1 



where-, 



(24) 



Note that the coefficients Bkn have to be calculated 
recursively since they contain coefficients B_, n _i. The same 
methodology can also be applied to find the square of dt to 



20 find: 



dt 2 - [ZAP 



2k - 2 \ = c x + c 2 



(p 2 ±B kl p 2 ™ 



p 4y LB k2 p 2k - 2 1+ 



{p't 



dx 1 



dx l 



B k3 p 2k ~ 2 | + « 



(25) 



23 



WO 2005/017563 



PCT/GB2004/003540 



where, 



(26) 



Hence, written out up to 6 th order in p, equation (25) 
becomes : 



A + AlP 2 + AP* + AP 6 + — = q + c 2 5 n p 2 + (c 2 5 21 + c 3J B 12 )/? 4 + (c 2 £ 31 + + 



10 



Thus the first two coefficients, using equations (18a-c),are 



r 



" i 

2>J 



Vs, 



Vp, 

1 



r.2 



0 » 



(27a) 



\ 



c„ = 



(27b) 



'11 



/ras 



&=1 
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Note that coefficient ci (equation 27a) can be directly 
interpreted . as the square of the difference in traveltime 
between the P- and PS-converted wave at vertical incidence 
(p=0 s/m) . The coefficient c 2 can not be interpreted so 
easily, although the numerator is equal to dto, but it has 
the same form as was found by Tessmer and Behle for PS- 
converted waves and hence it can be equaled to the inverse 

* 

■ 

of the square of an rms -velocity. 



25 



A Dix-Krey relation for reflection seismic receiver 
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Analogous to the derivation of a Dix-Krey type relation for 
the two-term traveltime difference approximation in the 
slowness domain, Dix-Krey type relations can also be derived 
for the approximations of the traveltime difference as a 
function of the difference in horizontal distance traveled 
by the P- and PS-wave. This is the topic of the current 
section. As will be seen, the result will relate the rms- 
velocities defined in the previous section to the products 
of P- and S-wave interval velocities. Starting from the 
squared rms-velocity : 



_ XLAfo* -*»*) 



vL = " — • (28) 



and proceeding by substituting equation (7). for the 
thickness o'f a layer k and the ratio of one-way vertical 
travel times (equation 8) into equation (28) results in: 

20 ELA(^-^*)=S*.,i(^* r * +v^Xvs k -vp k )=Y k M-<)-y s >yp>< • < 29 > 



Hence, combining equation (28) and (29), for layer n and n-1 
respectively gives : 

<vL,„ = Y k M W3r* > (30a) 
< _1 vL. n -i = T^te -<)-Vp k Vs k . (30b) 



Subtracting equation (30b) from equation (30a) leads to the 
product of P- and S-wave interval velocities in layer n: 

30 
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T/„ T/c u rms,n u rms,n—I 

Vp " s - — w^i — • <31) 



Now finally, using that: 



(r„"-<)=<-^- 1 , (32) 



equation (31) can be written: 



T/»-, T/ 0 . u rms,n o rms,n—\ 

Vp " Vs - (<-*;-) • <33) 

Equation (33) is the desired Dix-Krey type relation. It 
shows how the product of P- and S-wave interval velocities 
for a layer k, can be calculated once the rms velocities and 
the vertical incidence time differences for that layer and 
the previous layer are known. Equation (33) has exactly the 
same form as found previously for PS-converted waves by 
Tessmer and Behle (1988) . Note however that the vertical 
incidence traveltime differences and rms velocities are 
defined differently, as shown in the previous section. 



In the following it is illustrated how the derived 
approximations may be used with some numerical examples. 



1: Moveout correction using the two -term 
approximation in the slowness domain 

Synthetic data was calculated using a 2D reflectivity code 
for a horizontally layered model consisting of 6 homogeneous 
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layers overlying a half space. The medium properties are 
summarised in Table 1. In the "synthetic experiments, 801 
receivers were positioned symmetrically about the source 
position with a receiver-receiver spacing of 3.125m. Thus, 
.5 the maximum source-receiver offset considered here is 1250m. 
The sampling rate was 2ms. The source in the modeling was a 
P-wave point source, positioned on the surface, emitting a 
50Hz Ricker wavelet. As pre-processing . step the data were 
transformed to the time-slowness (i,p) domain. The pre- 
10 processed reflectivity data is shown in FIG. 3. 

■ 

The receiver functions, calculated by stabilised 
deconvolution, for this data are shown in FIG. 4A. The 
theoretical traveltime differences between the PS-converted 

15 waves and primary P-ref lections are calculated by ray- . 
tracing and are shown in blue. The two- term approximation 
-(eq. 5) -using exact pseudo rms-velocities is shown in cyan. 
It is clear that the two-term approximation is exact at zero 
slowness (p=0 s/m) . Moreover, the accuracy decreases with 

20 both increasing slowness and depth. This is because the 

* 

velocity increases with depth and the accuracy of the 
expansion is directly related to the product of the velocity 
and the slowness in . each layer. In FIG. 4B, the receiver 
functions have been NMO corrected using the two-term 
25 approximation and the linearly interpolated pseudo rms- 
velocities. It is clear that all identified events have been 
significantly flattened. The predicted residual moveout, 
calculated by extracting the two-term approximation from the 
ray- traced (exact) traveltime differences, is shown in blue. 

30 

In FIGs. 5, the area denoted by the black boxes in FIGs.4s 
is shown in detail. The two main events correspond to the 
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relative move out between the P-reflection and its PS- 
conversion at the fifth and sixth interface. Again, it is 
clear that the events have been significantly flattened 
using the two- term approximation. 

5 

Note that the moveout correction using the two- term 
approximation obviously requires a pseudo rms-velocity 
model. In the example here, the exact pseudo rms-velocities 
is calculated using eq. 6b and Table 1. The two-term 
10 approximation of the traveltime difference makes the problem 
amendable to velocity analysis (Taner and Koehler., 1969) by 

* 

reducing the number of parameters to two . 



Layer 


P-velocity [m/s] 


S-velocity [m/s] 


Depth [m] 


1 


1650 


400 

* 


0 


2 


1775 


700 


50 


3 


* 

1900 


950 


125 


4 


2000 


1000 . 


225 


5 


2250 


1125 


425 


6 


2750 


1375 


700 


7 


3000 


1500 


1300 



Table 1 Medium properties of the 6-layer over a halfspace 
15 model used in the reflectivity calculations. 

Example 2: Moveout approximations of the traveltime 

■ 

difference as a function of the difference in horizontal 
20 distance traveled 

In this example we use the exactly the same reflectivity 
data as was used in example 1, with the exception that the 
input data is not transformed to the time-slowness (x,p) 
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domain. Instead, as proposed in the section on series 
approximation of the traveltime difference as function of 
the difference in horizontal difference, -we calculate a new 

♦ 

type of receiver function formed by two-dimensional (2D) 
5 deconvolution of the vertical component (upgoing P-waves) 
out of the radial component (upgoing S-waves) data in the 
(x,t) domain. Two-dimensional deconvolution in the (x,t) 
domain is equivalent to division of each point in the 
corresponding f requency-wavenumber (f,k) ' domain. Hence, 
10 individual plane PS -wave components are shifted in time and 
space relative to the arrival time and location of the P- 
wave components with the same frequency and (horizontal) 
wavenumber . 

15 The result of this 2D-deconvolution is shown in FIG. 6A. The 
theoretical traveltime differences as function of the 
difference in horizontal travel distance, calculated by ray- 
tracing, are shown for each interface in blue. The two-term 
approximation, calculated from eq. (20) using coefficients 

20 eq. (27a, b) and exact medium properties, is shown in cyan. 
FIG. 6B shows (in green colour) successive higher-order 
approximations. The coefficients for the higher-order terms 
have been calculated from eq. (25), repeatedly making use of 
eq. (22), (24) and (26). Although it is not clear from FIGs. 

25 6 (or FIGs. 7), closer inspection shows that although the 
higher-order approximations are more accurate for small 
slownesses, they break down at a lower slowness. This has to 
do with the asymptotic nature of the series approximation. 

30. Although the 2D-receiver functions in FIG. 6A contain a lot 
of secondary events which cannot be explained by the ray- 
traced theoretical traveltime differences for each of the 

29 



WO 2005/017563 



PCT/GB2004/003540 



six interfaces, it is clear that those events that can be 
related to the interplay of a primary P-reflection and its 
PS-conversion, are approximated reasonably by the two-term 
approximation,, especially at small-slowness. 

■ 

In FIG . 7A, two events from the 2D-receiver function 
corresponding to the fifth and sixth interface, are shown in 
detail- FIG- 7B shows the corresponding zoom- in on the 
theoretical, two-term and higher-order traveltime difference 
curves - 

While the invention has been described in conjunction with 

# 

the exemplary embodiments described above, many equivalent 
modifications and variations will be apparent to those 
skilled in the art when given this disclosure. Accordingly, 
the exemplary embodiments of the invention set forth above 
are considered to be illustrative and not limiting. Various 
changes to the described embodiments may be made without 
departing from the spirit and scope of the invention. 



